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ABSTRACT 

Aims. Young stars interact vigorously with their surroundings, as evident from the highly rotationally excited CO (up to Eu/k = 4000 
K) and H2O emission (up to 600 K) detected by the Herschel Space Observatory in embedded low-mass protostars. Our aim is to 
construct a model that reproduces the observations quantitatively, to investigate the origin of the emission, and to use the lines as 
probes of the various heating mechanisms. 

Methods. The model consists of a spherical envelope with a power-law density structure and a bipolar outflow cavity. Three heating 
mechanisms are considered: passive heating by the protostellar luminosity, ultraviolet irradiation of the outflow cavity walls, and 
small-scale C-type shocks along the cavity walls. Most of the model parameters are constrained from independent observations; the 
two remaining free parameters considered here are the protostellar UV luminosity and the shock velocity. Line fluxes are calculated 
for CO and H2O and compared to Herschel data and complementary ground-based data for the protostars NGC1333 IRAS2A, HH 46 
and DK Cha. The three sources are selected to span a range of evolutionary phases (early Stage to late Stage I) and physical 
characteristics such as luminosity and envelope mass. 

Results. The bulk of the gas in the envelope, heated by the protostellar luminosity, accounts for 3-10% of the CO luminosity summed 
over all rotational lines up to 7 = 40-39; it is best probed by low-7 CO isotopologue lines such as C"*0 2-1 and 3-2. The UV- 
heated gas and the C-type shocks, probed by '"CO 10-9 and higher-/ lines, contribute 20-80% each. The model fits show a tentative 
evolutionary trend: the CO emission is dominated by shocks in the youngest source and by UV-heated gas in the oldest one. This 
trend is mainly driven by the lower envelope density in more evolved sources. The total H2O line luminosity in all cases is dominated 
by shocks (> 99%). The exact percentages for both species are uncertain by at least a factor of 2 due to uncertainties in the gas 
temperature as function of the incident UV flux. However, on a qualitative level and within the context of our model, both UV-heated 
gas and C-type shocks are needed to reproduce the emission in far-infrared rotational lines of CO and H2O. 

Key words, stars: formation - circumstellar matter - radiative transfer - astrochemistry - techniques: spectroscopic 

1. Introduction K), along with ro tationallv excited H?0 (up to 470 K) and 

OH ( up to 620 K; Ceccarel li etaP [19991 iGiannini et all Il999l 
In the embedded phases of low-mass star formation (Stages 206^ iNisini et al.li2002 ). The origin of this hot gas has been 
and I; Whi tney et al. | | 2003| ; | Robitaille etalj |2006), the mate- debated ever since. The shape of the CO 6-5 spectra - a nar- 
nal suiTounding the protostar is exposed to energetic phenoni- emission feature on top of a broader one - and the pres- 



ena such as shocks and ultraviolet radiation (iShu et al 



. ISpaans etalJll995t iBachiUer & Tafallalll999t lArce et al 



19871: 
2007h . 



ence of a strong, narrow CO 6-5 emission feature point at a 
mixture of quiescent and shocked gas, with the quiescent gas 
Observationally, this leads to much higher intensities in hnes of Hkelv heated bv the protostellar UV field along the cavity walls 
rotationally excited carbon monoxide, water and other species dSpaans et al.lll9fl . Spatially resolved CO 6-5 maps, showing 
than would be possible from gas that is only heated through col- extended narrow emission, provide additional support for this 
hsions with warm dust in the bulk of the envelope. The ubiq- scenario (van Kempen et al. 2009). However, it remains unclear 
uity of surplus rotationally excited emission was first established ^ quantitative level to what extent UV radiation and shocks 
through ground-based observations of the CO_y^ =_,6-5 ti;ansi- g^^^ contribute to heating up the gas and powering the emission. 



tion (upper- level energy Eg/k - 120 K: ISchuster et al 

iHogerheiide et al.|ll998h . The Infrared Space Observatory (ISO) The Herschel Space Observatory (iPilbratt et al.ll20l6l) offers 

subsequently detected CO up to the 21-20 line (Ei,/k = 1300 a new set of tools to tackle this question. Amongst its first results 

are detections of highly rotationally excited CO (up to 7 = 38- 

* Herschel is an ESA space observatory with science instruments 37, E^/k = 4100 K) and H2O (up to 640 K) in the embed- 

provided by European-led Principal Investigator consortia and with im- ded low -mass protostars HH 46 and DK Cha dvan Kempen et al.1 

portant participation from NASA. l2010allbl) . A quantitative radiative transfer model of HH 46 was 
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successful at explaining the CO line fluxes with a combination of 
UV-heated gas in the outflow cavity walls and small-scale C-type 
shocks along the walls (van Kempen et al. 2010b). However, the 
model was unable to disentangle the origin of the observed H2O 
lines in a similar fashion. The goals of the cu rrent paper are (1) 
to present our model from Ivan Kempen et al.l in more detail; (2) 
to check how robust the model fit is for CO; (3) to derive a quan- 
titative fit for H2O; and (4) to extend the analysis to DK Cha and 
a third source (NGC1333 IR AS2A) for which CO and H2O data 
have since become available dKristensen et al.l2010HYildiz et al.l 
The three test cases are all nearby protostars (d - 180- 
450 pc) and cover the evolutionary range from deeply embedded 
Stage to late Stage I. 

The question of the origin of the hot gas observed with 
Herschel is not limited to embedded low-mass young stellar ob- 
jects (YSOs). Herschel observations of rotation ally excited CO 
and/or H2 O include intermediate-mass YSOs (Fichet al. 2010; 
iJohnstone et al. 2010). high-mass YSOs (.Chavarrfa et al..,2010,) . 



Orion Bar photon-dominated region dHabart et al.l 



circumstellar disks (iMeeus et aLll2O10 : Sturm et al. 



2010I) . the 
2OIO I), and 



even extraga lactic sources like the ultraluminous infrared galaxy 
Mrk 231 (Ivan der Werf et al]l201 0). A quantitative explanation 
of the origin of the observed emission is still missing in most of 
these cases. Our model may help in constraining possible excita- 
tion mechanisms in environments other than low-mass YSOs. In 
particular, the treatment of how the gas is excited by UV fields 
and shocks can be adapted for different types of sources. 

This work is part of the two complementary Herschel 
Key Programmes WISH and DIGIT. The primary goal of 
WISH, or "W ater in star-forming regions with Herschel" 
(Ivan Dishoeck et al. 2011 ), is to use H2O as a probe of the phys- 
ical and chemical conditions during star formation and to fol- 
low its abundance as a function of evolutionary phase. DIGIT, 
or "Dust, ice, and gas in time" (PI: N. Evans), aims to study the 
evolution of the gas during different stages of star formation, as 
well as that of the dust grains and their icy mantles. Observations 
of CO are an integral part of both programmes, because its sim- 
ple chemistry and easy use in radiative transfer codes make CO 
an excellent tool to constrain physical conditions and processes. 
Indeed, the approach in the current paper is to fit the free param- 
eters in the model to the CO data and then to assess how well 
the model reproduces the H2O data. More generally, the model 
is ultimately intended as a tool in using the Herschel data as a 
probe of the energetics and dynamics of star formation, in partic- 
ular how and where the young star deposits energy back into the 
parent molecular cloud. Future work in WISH and DIGIT will 
explore these questions in more detail. 

A detailed description of the model appears in Sect. |2] The 
Herschel data and complementary ground-based data are sum- 
marised in Sect. [3] The model results for CO and H2O are pre- 
sented in Sect.|4]and discussed in Sect.|5] Conclusions are drawn 
in Sect.lH 

2. Model description 

The far-infrared (far-IR) and sub-millimetre (sub-mm) molecu- 
lar line emission from embedded low-mass protostars has tradi- 
tionally been interpreted with purely spherical envelope mod- 
els, but such models cannot reproduce the observed intensi- 
ties of rotationally exci ted lines like CO 6-5 and higher (e.g., 
Ivan Kempen et alJl2009 ). Our model goes a step further by in- 
troducing a bipolar outflow cavity with two additional heating 
mechanisms: UV irradiation of the gas in the cavity walls and 
small-scale C-type shocks along the walls (Fig.[T]i. Rotationally 




Fig. 1. Schematic view of the physical components in an embed- 
ded Stage or I young stellar object: a passively heated envelope 
(yellow, with a hot core shaded red), a bipolar jet (green, with 
bow shocks and internal working surfaces), UV-heated outflow 
cavity walls (red), and small-scale shocks along the cavity walls 
(blue). Not visible on this scale is an embedded circumstellar 
disk with a radius of ~ 100 AU. 



excited emission may also arise from the bipolar jet and the asso- 
ciated internal J-type shocks, from gas heated by X-rays, or from 
a circumstellar disk embedded within the envelope. However, 
our goal is not to create a single all-inclusive model of an embed- 
ded protostar. Rather, the main question in this work is what the 
Herschel observations can teach us about the passively heated 
envelope, the UV-heated cavity walls and the small-scale shocks. 
Investigating the quantitative contributions from other compo- 
nents will be left for future studies. Furthermore, the current fo- 
cus is on the integrated line intensities only. With regards to the 
Herschel-PACS data, we are only interested in the emission from 
the central spaxel of the full 5x5 spaxel array. Future updates 
of the model will explore the spatial extent of the emission and 
could add dynamical processes to allow for the computation of 
line profiles. 



2.1. Passively heated envelope 

Low-y lines from CO and its isotopologues and similar low- 
excitation lines from other species are commonly modelled with 
a sph erical envelope (e.g.. |j0rgensen et aP 120021 : ISchoier et al.l 
12002 ). A spherical envelope is also the starting point of our 
model, since the warm inner region can contribute to emission 
in the higher-7 lines. The gas in the envelope is well coupled to 
the dust and is only heated passively or indirectly, i.e., the dust 
is heated by the protostellar luminosity and the gas attains the 
same temperature through collisions with the dust. 

The density and temperature profiles of each source are con- 
str ained from cont i nuum observation s. Following t he method 
of lJ0rgensen et al.l (|2002|) and Schoier et alj (l2002i) . the con- 
tinuum emission for a grid of envelope models is calculated 
with the ID continuum radiative transfer program DUSTY 
(Ivezic et al. 1999). The best-fit model is determined by com- 
paring the emission to spectral energy distributions (SEDs) and 
sub-mm brightness profil es compiled by FroebrichI (l2005 l) and 
iDi Francesco et al. I (l2008l) . The inner radius of the envelope is 
set at a dust temperature of 250 K (about 30 AU or 0.15" at 
200 pc; Table |2]l. If any material is present on smaller scales, 
its contribution to the overall line emission would be too diluted 
in the ~10" Herschel beam to have a significant impact on the 
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model results. The envelopes are optically thick at UV and vis- 
ible wavelengths, so most of the stellar radiation is reprocessed 
by the dust. For a constant stellar luminosity, the assumed stel- 
lar temperature therefore does not have a significant effect on 
the resulting dust temperatures or the molecular line emission 
(lJ0rgensen et al. 2002). This was verified by running part of the 
DUSTY grid for stellar temperatures of both 5 000 and 10000 
K. 

The remaining three free parameters are fitted in a pro- 
cedure: the power-law index for the density profile (n oc r '^), 
the optical depth at 100 fim (tioo), and the size of the envelope 
(Y = rout/nn)- For our line radiative transfer models, the en- 
velopes are terminated at the point where the density reaches 10 ^ 
cm or the temperature reaches 10 K (|j0rgensen et al.ll2002h . 
This point is denoted renv The best-fit values for the three test 
cases are listed in Table |2] in Sect. [3] and the associated density 
and temperature profiles are plotted in Fig. |5] 

The Doppler broadening parameter b for the passive enve- 
lope models is fixed at 0.8 km s"\ a typical value for the qui- 
escent gas in low-mass YSOs as determin ed from optically thin 
C"*0 observations (iJ0rgensen et al.l2002l) . The ongoing collapse 
of the envelope onto the star and di sk is s imulated with a power- 
law infall velocity profile v oc r^^/^ ( IShul[T977l) . scaled to 2.0 km 
s"' at the inner edge of the envelope. Details on how the abun- 
dances are obtained and how the molecular line intensities are 
calculated are provided in Sects. I272!2l and l272. 31 



2.2. Photon-heated cavity walls 
2.2.1 . Density and temperature 

The second model component is the UV -heated gas in the 
walls of the outflow cav ities (red in Fig. [Ij ISpaans et~anil995t 



Ivan Kempen et al. 200^. For an ambient cavity density of a few 
10^ cm (Neufeldet al] |2009l) . the total gas column from the 
central source to the point where the cavity intersects the outer 
edge of the envelope (a distance of ~ 10'^ cm) is a few 10^" cm"^. 
This corresponds to a UV optical depth of a few tenths, so UV 
photons produced near the protostar can freely reach the cav- 
ity walls. When they do, they produce a photon-dominated zone 
where the gas is heated to temperatures of a few hundred K. 

The first step in modelling the emission from the UV-heated 
gas is to carve out a bipolar cavity from the spherical env elope, 
similar to the procedure followed by Bruderer et al.l (l2010l) . Mid- 
infrared Spitzer Space Telescope images of HH 46 and similar 
YSOs like L 1527 show ellipsoidal cavities, with a ratio between 
the semimaior and semiminor axes (ocav and bc^w) of about 4 
dWusamy et al. 2007; Tobin et al. 2008). Large-scale CO 3-2 
maps of NGC1333 IRAS2A, t racing the swept-up outflow gas, 
show roughly the same s hape (ISandell et £0119941) . DK Cha is 
oriented close to pole-on (iKnedl 19921) . so the shape and size of 
its outflow are unknown. 

Tests of our model show that the variation in the computed 
CO and H2O line fluxes resulting from changes in the shape and 
size of the cavities is within the overall model uncertainties, con- 
sistent with the more detailed tests of Bruderer et al. ( 2010). Our 
approach is therefore to fix the ratio between the outflow axes 
(acav/^c av) of all thrcc sources to the value of 4.3 measured for 
HH 46 (iVelusamv et alj|2007l) . and to scale the absolute length 
of the axes with the envelope radius. In other words, the ratios 
dcav/fenw and fecav/'"env of all three sources are fixed to the val- 
ues of 2.1 and 0.50 derived for HH 46. The envelope radii and 
outflow axis lengths are summarised in Table |2] 




Fig. 2. Schematic view of the envelope and the bipolar outflow 
cavity. In each quadrant, a point on the cavity wall is uniquely 
identified by its distance rcav to the central star at the origin (Eq. 
©). 



The two cavities for each source are aligned along the z axis, 
with the ends meeting at z = (Fig. |2]l. The ellipsoidal cavity 
wall in the upper right quadrant of a cylindrical coordinate sys- 
tem is defined as 




(1) 



The points along the wall can also be characterised by their dis- 
tance Tcav to the protostar. 



+ z^ 



(2) 



which is used for some figures in Sect. |4] 

The only UV source included in our model is t he protostar. 
UV photons can also be pro duced at the bow shock dBohm et al.l 
ll993l:lRavmond et al.' 1997) or by small -scale shocks within the 
cavity or along the walls (Curiel et al. 1995; Saucedo et al. 20031 
Wa lter et"al]|2003h . but these contributions are highly uncertain. 
The protostellar UV luminosity itself is also highly uncertain, 
and is therefore treated as a free parameter Classical T Tauri 
stars have typical L^y of 10~^ '^-10~''^ times the total stellar lu- 
minosity (llnglebv et al.ll20lTl) . Scaling up to the higher accre- 
tion rates found in embedded YSOs, anything from 0.01 to 1 Lq 
would be a reasonable value for our model. The best-fit UV lu- 
minosities derived in Sect. l4.2l range from <0.03 to 0.08 Lq. An 
Luv of 0.1 Lq would give an unattenuated UV flux at 100 AU 
of Go ^ lO'* relative to the mean interstellar radiation field of 
lHabingl (ir96{^. 

Due to the UV field impinging on the cavity walls, the gas 
is heated to higher temperatures than it attains in the passively 
heated, purely spherical envelope. The temperature can be cal- 
culated by solving the heating and cooling balance as function 
of the physical conditions. Many such models are available in 
the literature, but they show substantially different results, some- 
times even at the qualitative level (Rollig et al. 2007). The con- 
sequences of these uncertainties for our model results are dis- 
cussed extensively in Sect. 15.11 For the initial results in Sect. |4] 



3 



R. Visser et al.: Modelling Herschel observations of hot molecular gas from embedded protostars 



the PDR surface temperature grid of iKaufman et alJ (Il999h is 
adopted as our standard prescription. This grid was calculated 
for densities up to 10^ cm"^; for higher densities, it is extrap- 
olated as 1/«(H2) based on the expected dominant heating and 
cooling mechanisms (Sect. 15. 11 1. 

The iKaufman et al.l (1 19991) grid as such can only be used 
at the cavity wall. At any depth into the envelope, the tem- 
perature needs_to_be corre cted for attenuation of the UV field 
by dust. iRoUig et"al] (l2007b plotted the temperature as function 
of the visual extinction. Ay, for several densities and unattenu- 
ated UV fluxes. In each case, the temperature decreases roughly 
as exp(-0.6Av). In order to use this depth dependence in our 
model, the amount of extinction towards any point (R, z) is cal- 
culated with the continuum r adiative transfer code RADMC 
dPuUemond & Dominikll2004h . The input stellar spectrum is a 
10000 K blackbody scaled to a given UV luminosity between 6 
and 13.6 eV. This input spectrum is a reasonable approximation 
of the real spectrum emitted by a low-mass protostar w ith a UV 
excess due to accretion (e.g.. Ivan Zadelhoffet ani2003h . The lo- 
cal UV flux from RADMC, f radmc (also integrated from 6 to 
13.6 eV), can be considered the product of the unattenuated flux 
and an extinction term, 

^ unatt exp(-Tuv) , (3) 

where Funatt is the UV flux through a unit surface arising from 
a source at a distance V^^+z^ xhe UV optical depth, tuv, is 
effectively an average over the many possible paths a photon can 
follow from the source to the point ( R,z). It is converte d into a 
visual extinction via Ay = ruv/3.02 (Boh lin et al.lll978h . Along 
with the UV field, RADMC also calculates the dust temperature 
along the cavity wall and throughout the envelope. This is used 
instead of the dust temperature from the ID DUSTY models. 

One caveat with RADMC is that it treats scattering of the 
radiation by dust as an isotropic process, while in reality grains 
are partially forward scattering (Li & Draine 2001). Tests with 
the radiative transfer code of van Zadelhoff et al. (2003), which 
does allow for anisotropic scattering, show flux differences of at 
most 20% between an isotropic scattering function and the scat- 
tering function of iLi & Drainei This is within our overall model 
uncertainty. 

The present goal is to reproduce the observed integrated line 
fluxes, not the detailed line profiles. Hence, the Doppler broad- 
ening and infall velocities for the UV-heated component are kept 
the same as in the passively heated component. 

2.2.2. Abundances 

Several thousand grid points are required to properly sample the 
envelope and the cavity walls in the 3D radiative transfer cal- 
culations (Sect. 12. 2.31 1. It would be too time consuming to cal- 
culate the abundances of CO and H2O for each individual point 
with a traditional chemical code. Inst ead, the ab undances are ob- 
tained according to the method of Bruderer et al. (2009). A grid 
of chemical models is pre-computed for the appropriate range 
of densities, temperatures, unattenuated UV fluxes and extinc- 
tions. Abundances for each individual point in the envelope or 
the cavity wall are interpolated from this grid when constructing 
the source models for the radiative transfer code. 

The b asis of our chemical r eaction network i s the UMIST06 
database (Woodal l et al.l l2007h as modified by iBruderer et alJ 
(2009), except that X-ray chemistry is not included. The network 
allows all neutrals to freeze out onto the dust to account for the 
depletion of many gas-phase species in the colder parts of the 



Table 1. Initial abundances of molecules in our chemical grid. 



Species 


Abundance 


H2 


1.00(0) 


CO 


1.57(-4) 


HoOice 


2.03(-4) 




4.94(-5) 



Notes. The notation a(b) means ax 10*. Abundances are given rel- 
ative to H2. Initial abundances of elements present in the UMI ST06 
database but not listed in this table are as in lAikawa et al.l ( l2008h and 
.Bruderer etal.,aOO^) . 



envelope. Molecules can return to the gas phase through ther- 
mal desorption and photodesorption. Binding energi es and pho- 
todesorption yields for CO and H2O are taken from Fraser et 
(2001), Bisschop etal] ( |2006|) and lOberg et"a l. (2007, 20oC 
Photodesorption occurs even in regions of high extinction be- 
cause of a cosmic-ray-induced UV flux of about 10'* photons 
cm~^ s~', equiva lent to 10 times the mean interstellar flux 
(IShen et a l.''2004'). In the gas phase, UV photons can dissociate 
or ion ise m any species; the relevant rates are calculated accord- 
ing to lvanDishoecket al.1 ( 120061) for a 10 000 K blackbody stel- 
lar spectrum, appropriate for a low-mass protostar with a UV ex- 
cess. Self-shielding of H2 and CO, which causes their dissocia- 
tion rates to decrease more rapidly than those of other spe cies, is 
included according to Draine & Bertoldi ( 1996) and Visse r et aH 
(I2009|)^ The cosmic -ray ionisation rate of H2 is set to 5 x 10 
s-^ (lDalgarn6ll2006h . 

The grid of pre-computed abundances, from which abun- 
dances for the three sources are interpolated, is obtained by ad- 
vancing the chemical reactions for a period of lO*" yr. As starting 
conditions, all carbon is locked up in CO, the remaining oxygen 
in H2O, all nitrogen in N2, and all remaining hydrogen in H2. 
Reflecting the most likely formation pathways, CO and N2 start 
in the gas phase, while H2O starts frozen out on the dust. H2 is 
predominantly formed on the grains, but because it rapidly evap- 
orates even at 10 K, it starts in the gas phas e. E lemental abun- 
dances are adopted from Ai kawa et al.l (12008 *) and'Bruder er et al.l 
(2009); the initial composition is summarised in Table[T] 

2.2.3. Radiative transfer 

Level populations and line intensities for the passively heated 
envelope and the UV-heated outflow cavity walls are cal- 
culated with the new full- 3D radiative transfer code LIME 
(iBrinch & Hog erheiide 2010). LIME uses an irregular grid with 
points sampled randomly but weighted according to a user- 
defined criterion such as density or temperature. This naturally 
results in a finely spaced grid in the regions most important for 
the astrophysical problem at hand. 

The source geometry of our models introduces a substan- 
tial difficulty: the high-excitation CO and H2O lines originate 
in a thin layer of hot gas along the cavity wall, so an accurate 
flux calculation requires a high grid point density inside and im- 
mediately outside this small region. Furthermore, the grid has 
to be able to resolve emission on scales ranging from ~10 to 
~10^ AU. This cannot be accomplished with the default density 
or temperature weighting functions. Extensive convergence tests 
were caiTied out to find an alternative grid sampling scheme. 
The optimised grid consists of 30 000 points, picked randomly 
in logr between and r^nw For every point with radial coor- 
dinate r, a point (^cav,Zcav) is located on the cavity wall such 
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that V-'^cav + -^cav - '"■ Let ycav be the angle between the point 
(^cav,Zcav) and the midplane, i.e., tanycav - Zcav/^cav The angu- 
lar coordinate y of each point is then determined as follows: 

1 . one third of the points are placed exactly at the cavity wall 

(r = Tcav); 

2. one third of the points are picked with a random angle y such 

that 0.9 < r/rcav < 1; 

3. one third of the points are picked with a random angle y such 
that < y/ycav ^ 1 and then weighted by ^/yJy^, which 
favors points near the midplane over points near the cavity 
wall. 

The grid thus co nstructed is smoothed by five iterations of 
Lloyd's algorithm (lLlovdlll982h . which moves each point some- 
what towards the centre of mass of its Voronoi cell. The purpose 
of this smoothing step is two-fold: (1) to turn some "needle- 
shaped'' cells_created by the random setup into more regular 
shapes dBrinch & H ogerhei id3l2010h : and (2) to push some grid 
points from the first sampling group into the cavity wall, so that 
the edge between the envelope and the cavity is well sampled on 
either side. The first two groups of points together cover the re- 
gion of hot gas responsible for the more excited lines. The third 
group of points covers the bulk of the envelope, where the low-7 
emission originates. 

With the grid established, LIME employs the standard two- 
step method for problems where local thermodynamic equilib- 
rium (LTE) does not apply. The first step is an iterative proce- 
dure to calculate the rotational level populations based on the 
balance between radiative excitation, collisional excitation and 
collisional de-excitation. The second step is the ray tracing, i.e., 
calculating the emergent spectrum at a given distance and incli- 
nation. The collisional rate coefficients required for the first step 
are collect ed from the Le iden Atomic and Molecular Databasqj 
(LAMDA: ISchoier et al .''20Q5). where the most recent d ata files 
for CO and H2O are from Yang et al. (2010) and Fa ure et alJ 
(l2007h . respectively. The ortho-to-para ratio for H2 is taken to 
be thermalised, with a maximum value of 3. The ortho-to-para 
ratio of H2O is fixed at 3. Dust is included in the model with 
a standard gas-to-dust ratio of 100 and the OH5 opacities from 
[Ossenkopf & Henning. (1994), appropriate for dust grains with a 
thin ice mantle. 

Images are created at (X'l resolution (18^5 AU at the dis- 
tances of the three sources; Table |2]l to get a sufficiently fine 
sampling of the hot inner regions. In order to compare the 
model results to the observations, each image is convolved 
with a Gaussian beam profile. Beam sizes are calculated as the 
diffraction-limited beam for the wavelength and the telescope 
(Herschel, APEX or JCMT) at which the simulated line was ob- 
served. In the case of lines observed with Herschel-PACS, fluxes 
are extracted from the central 9'.'4x9'.'4 to mimic the instrument's 
central spaxel. For a given density, temperature and abundance 
structure, the estimated uncertainty on the CO line fluxes ranges 
from 30% for the 2-1 line to a factor of 2 for the high-7 end of 
the ladder For H2O, it is a factor of 2 for all lines. The uncer- 
tainty is mainly due to the complex geometry and the small size 
of the emitting region for the high-7 lines. The source models 
themselves lead to uncertainties of similar magnitude in the line 
fluxes, mostly because of the problems with the gas temperature 
discussed in Sect. 15. II 

Integrated intensities from APEX, the JCMT and Herschel- 
HIFI ( J rmbdu in K km s ' ) and integrated fluxes from Herschel- 



PACS ( / Fydv in W m-2) can be mutually converted through 



r ikD. r 

I Fydv = -p- I Tmbdy , 



(4) 



with Q(/l) the solid angle subtended by the beam at the given 
wavelength. Denoting the full width at half maximum (FWHM) 
of the beam as ©(A), the sohd angle for APEX, flie JCMT and 
HIFI is (n/A In 2)©^. The formula is more complicated for PACS, 
because the size of the spaxels has to be taken into account. For 
just the central spaxel, the solid angle is 



f^PACsC'^) 



9'.'4 



l2 



erf(9'.'4 y/\n2/@) 



(5) 



with erf the error function and @{A) the FWHM of the 
diffl-action-limited Herschel beam. 



2.3. Small-scale shocks 

The third model component is a series of C-type shocks mov- 
ing along the cavity walls. These shocks may be powered by 
the wind launched from close to the protostar or th e surface of 



the inner disk, or by the expansion of the iet (e.g..lLadalll985 



Shang et all 120061: lArce & Sargen t 2006; Santiago- Garcfa et al.1 



http://www.strw.leidenuniv.nl/~moldata 



20091) . Shocks are created when the wind hits the outflow cavity 
walls, heating the gas to temperatures of more t han 1000 K. 

Th e model procedure is the same as that of iKristensen et al] 
(l2008h . except for a change in the geometry. In the original pro- 
cedure, ID shock model results were pasted along a parabola 
to simulate emission from a bow shock observed in H2 in the 
Orion Molecular Cloud. Here, in order to simulate C-type shocks 
along the cavity walls, the ID shocks are pasted along the el- 
lipsoids from Sect. 12.2.11 The shock velocity is assumed to be 
constant along the entire length of the wall and the pre-shock 
density is that of the envelope at the given distance from the 
protostar The magnetic field is set to yS[2n(H2)/cm"-']'^^, with 
/5 fixed at 1 /iGauss. Abundances as function of depth into the 
shock are calc ulated from a lim ited chemical network centred 
on H2O (Wa gner & Grafiiri987h . which does not include pho- 
todissociation or other UV-driven reactions^ 

Th e shock model results of iKaufman & Neufeldl 
(1996) a re used in co mbination with the model results of 
Flower & Pineau des F orets (2003) and Kristensen et al. (in 
prep.) to compute the emergent line intensities. The first step is 
to estimate the extent of the line-emitting region by measuring 
the FWHM of the cooling-rate profile of the relevant molecule 
through the shock. The bottom panel of Fig.[3]shows the cooling 
profiles from Flower & Pineau d es Forets for CO and H2O, for 
a shock velocity of 20 km s"' and pre-shock densities of lO'* 
and 10^^ cm"^. The corresponding temperature and abundance 
profiles are plotted in the top panel. CO and H2O cool the 
same part of the shock and their cooling lengths are nearly the 
same. Furthermore, the cooling lengths are nearly inversely 
proportional to the pre-shock density, i.e., the solid and dotted 
profiles look similar but are shifted on the logarithmic scale. The 
cooling lengths provide an estimate of the extent over which the 
shock is emitting, and as such are an effective measure of the 
beam-filling factor 

The second step is to divide the cavity wall into discrete 
segments according to the envelope density. The upper limit of 
the Kaufman & Neufeld ( 1996) grid is 10*^ cm"^; the pre-shock 
density for the denser parts of the envelope is set to that value. 
Figure |4] shows the adopted segmentation for HH 46; the width 
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Fig. 3. Top: Temperature of the neutral gas (black) and abun- 
dances of CO (red) and H2O (blue) as function of depth into 
a shock for a shock velocity of 20 km s ' and pre-shock den- 
sities of 3 X 10'' (soHd) and 1 x 10"* cm"^ (dotted). Bottom: 
Normalised cooling rates of CO and H2O for the same two 
sets of conditions. All data are taken from the models of 
iFlower & Pineau des Forets (.2003.) . 



of the coloured area is a rough representation of the size of the 
emitting region. The CO an d H2O line intensities are estimated 
from the results tabulated bv lKaufman & Neufeldl The emission 
is assumed to fill each of the segments, thus creating an emis- 
sion map for each line. Just like the raw images produced by 
LIME (Sect. 12.2.31 ). the maps are convolved with a Gaussian 
beam profile of the appropriate diameter For lines observed with 
Herschel-PACS, the flux is extracted from the central 9'.'4 x 9'.'4. 
The fluxes are not corrected for extinction by dust in the enve- 
lope. 

Our method is specifically designed to simulate optically thin 
emission from shocks. Low- J CO lines typically observed in and 
associated with outflows, such as the 2-1 and 3-2, originate in 
cold (~100 K) swept-up gas present on larger spatial sca les than 
the shocks along the cavity wall (e.g.. lArce et al.ll2007h . These 
low-7 lines are often optically thick and only probe the surface 
layer of the outflow, where the velocity gradient is largest. Since 
our focus is on the hot gas observed with Herschel, no attempt 
is made to quantify the emission from the colder swept-up gas. 

An important caveat in the shock models is that they were 
calculated without a UV field. Several groups are working on 
irradiated shock models featuring a self-consistent treatment 
of the UV field in both the pre- and the post-shock gas (e.g., 
lLesaffTeetal.li201li: M. J. Kaufman, priv. comm.), but no such 
models are currently available. Qualitatively, a UV-irradiated 
shock is expected to be smaller and hotter, and to have a lower 
H2O abundance. The latter is a direct effect of the photodisso- 
ciation of H2O. The smaller spatial extent and the higher peak 
temperature are due to the higher degree of ionization in the irra- 
diated pre-shock gas: the ion-neutral coupling length decreases. 
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Fig. 4. Segmentation of the outflow cavity wall of HH 46 into 
regions of different pre-shock densities. The cooling length is 
larger for lower densities, so the shocks are wider at larger radii. 



so the shock front is compressed and the same amount of en- 
ergy is imparted over a shorter distance. The shock results are 
discussed in Sect. |4] with these qualitative arguments in mind. 

3. Source sample and observations 

The three test cases for the model are the low-mass 
YSOs NGC1333 IRAS2A, HH 46 IRS and DK Cha (IRAS 
12496-7650). Each source has been studied extensively in 
the past, inc luding recent observation s with PACS and HIFI 
on Herschel ("van Kempen et al]|20I0alfbl : lKristensen et al.ll2010l : 
[Vildiz et al. 2010). The availability of Herschel data and a large 
body of complementary data is one reason to choose these three 
sources. Another is that they offer various challenges to the 
model by spanning a range of evolutionary stages and orienta- 
tions: IRAS2A is a deeply embedded Stage source seen nearly 
edge-on (Sandell et al. 1 994), HH 46 is a Stage I source seen at 
a 53° incHnation (Nishi kawa et al.ll2008l) . and DK C ha is in tran - 
sition from Stage I to II and is seen nearly pole-on ( Kneel[T992h . 
The basic observational characteristics and model parameters for 
each source are summarised in Table|2l 

Far-IR spectra of t he three sources were obtained with PACS 
and HIFI on Herschel dPilbratt et alJl2010tlPoglitsch et alj|2010l: 
Ide Graauw et a l. "20T0I). The HIFI data used in this work were 
pubhshed bv [Kristens en et all 12010) and lYildiz eTaP (l2010h . 
The PACS spectra of H H 46 and DK Cha were published by 
Ivan Kempen et al.1 (l2010a.b) : the PACS spectrum of IRAS2A is 
previously unpublished. All PAC S data wer e rereduced with the 
standard pipehne in HIPE v6.1 ( IOttll2010h . A spectral flatfield 
was applied to the data to improve the signal-to-noise ratio. 

PACS consists of a 5x5 integral field unit with 9'.'4 x 9'.'4 
spaxels. For HH 46, the fluxes of the central 9'.'4 region are 
measured directly from the central spaxel and are accurate to 
10-20%. The observations of DK Cha and IRAS2A were mis- 
pointed by a few arcsec. For both sources, the line emission is 
concentrated near the location of the continuum emission. In or- 
der to measure the line fluxes, the spectra in the four spaxels 
with the brightest continuum emission were first summed. The 
ratio of the continuum emission in these four spaxels to the to- 
tal continuum flux in the 5x5 array was then used to obtain a 
wavelength-dependent curve to convert the four-spaxel extrac- 
tion into a flux in the entire field of view. The fluxes were mul- 
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Table 2. Observational characteristics, model parameters and derived properties for the three sources. 



Source 


d 




Y 


P 


T 100 


r-m 


env 




n(''cnv) 




^cav 


^cav 


i 








(pc) 










(AU) 


(AU) 


(cm-3) 


(cm-3) 


(Mo) 


(AU) 


(AU) 


n 


(Lo) 


(km s-') 


IRAS2A 


235 


37 


500 


1.5 


0.4 


28.0 


14000 


1.5(8) 


1.3(4) 


2.0 


29 000 


7 000 


90 


<0.03 


15 


HH46 


450 


26 


3000 


2.0 


1.9 


34.6 


16100 


1.1(9) 


5.0(3) 


1.7 


34000 


8 000 


53 


0.08 


15 


DK Cha 


178 


36 


3000 


2.1 


0.1 


25.6 


5 700 


8.5(7) 


1.0(3) 


0.02 


12000 


2 800 





0.05 


>20 



Notes . References: Kristensen et al. (subm.) for L^,: ISandell et al. l 11994) and 'Hirota et al.' ("2008^ for N GC1333 IRAS2A; iGraham & Heved 
( ll989h , IVelusamvetal.l ( 120071) and lNishikawa et alj l l2008l) for HH46; Knee ( 1992) and Whittet et al] ( 119971) for DK Cha. 
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Table 3. Measured line fluxes (10 W m ^) from the central 
region around each source. 
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Fig. 5. Density and temperature profiles for the passively heated 
envelopes of the three sources. 



tiplied by the fraction of the flux that falls in the central spaxel 
(ranging from 70% below 100 /im to 40% at 200 jum) to esti- 
mate the total flux that would have been in the central spaxel 
of a well-pointed observation. These corrections assume a point 
source and may therefore overestimate the flux in the central 
spaxel. Excluding this uncertainty, the fluxes for IRAS2A and 
DK Cha are accurate to about 30%. Fluxes and upper limits for 
all three sources are tabulated in Table[3] 

The Herschel data are combined with observations of 
lower-7 CO and CO isotopologue lines (up to 7-6) obtained 
with the James Clerk Maxwell Telescope (JCMT) and the 
Atacama Pathfind er Experiment (APEX; J0rgens en et al. 200'5 
Ivan Kempenetali 2006. 2009; Yildiz et alBoH) . The diff^erent 
beam sizes for the JCMT, APEX and Herschel are taken into 
account by convolving the line fluxes from the model with a 
wavelength-dependent beam size appropriate for the instrument 
with which each line was observed. 



4. Results 

4.1. Physical and cinemical structure 

The best-fit DUSTY envelope model of HH 46 has a radius of 
16 100 AU (cut off' at 10 K) and a mass of 1.7 Mq (Table EJ. 
The envelope of NGC1333 IRAS2A is a httle smaller (14 100 
AU), but because of the shallower density profile its mass comes 
out 20% higher at 2.0 Mq. DK Cha has the smallest and least 
massive envelope: 5700 AU and 0.02 Mq. The density and tem- 
perature profiles of each envelope are plotted in Fig.|5] 

When the cavities are carved out of the spherical envelopes, 
the gas and dust are heated up by the stellar UV flux. The best- 
fit UV luminosities for the three central sources are determined 
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28-27 


2240 


93.35 
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74.89 


<29 
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34 
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Notes. Typical uncertainties are 30% for IRAS2A and DK Cha and 10- 
20% for HH 46. Upper limits are given at the 3cr level. Lines in italics 
are blended (CO 31-30 with an OH line) and their fluxes are treated as 
upper limits. Since the observations are compared with model images 
that are convolved with the Herschel point-spread function, the fluxes 
in this table are non-PSF-corrected. 
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Fig. 7. Gas temperature profiles due to UV heating of the outflow cavity walls. Shocks are not included. The bottom row zooms in 
on the regions indicated in the top row. 

tances, the decreasing density allows the gas to become hotter 
than the dust, reaching a maximum of 260 K at 670 AU. 

The full 2D gas temperature profiles in Fig. Q show that the 
effects of UV heating are not limited to just the cavity wall, 
but penetrate some depth into the envelope according to the 
exp(-0.6Av) depth dependence from Sect. l2.2.TI The thickness 
of the UV-heated layer ranges from about 100 AU at r = 1000 
AU to about 2000 AU at the outer edge of the envelope. Hence, 
the UV-heated layer is about a factor of 10 thicker than the 
shock-heated layer (Fig.O. Accounting for the predicted effects 
of the UV field on the shock structure, the thickness ratio be- 
tween the two layers would be even larger (Sect. l273] ). 

In the purely spherical models for IRAS2A and HH 46, the 
CO abundance follows a drop profile (top panel of Fig. [8] and 
top left panel of Fig. |9]): it is high (~10 in the warm inner 
parts of the envelope, low (~10"*-10"^) at larger radii due to 
freeze-out onto cold dust, and high again (~10""*) towards the 
outer edge because of the low density and corresponding long 
freeze-out timesca le (ISchoier et alJl2004t lJ0rgensen et alJl2005E 
lYildiz et al. DK Cha has a somewhat warmer envelope 

(Fig. |5]l, preventing CO from freezing out and keeping its abun- 
dance at 10 at all radii. 

There is no freeze-out either along the cavity wall in any of 
the 2D models (bottom panel of Fig. [8] and bottom left panel of 
Fig.|9|l, because the dust temperature now exceeds the CO evapo- 
ration temperature everywhere (Fig.|6]l. The 10 000 K blackbody 
adopted for the stellar spectra is energetic enough to dissociate 
some CO in the cavity wall, reducing its abundance to 10"''-10"^ 
in IRAS2A and HH 46, and -10"^ in DK Cha (Fig.E bottom). 
Moving from the cavity wall towards the midplane, the increas- 
ing amount of extinction shields CO from the dissociating pho- 
tons, and the abundance goes back up to ^lO"'* (Fig.|9j see Fig. 
[To] in the online appendix for the other two sources). The dust 
near the midplanes of HH 46 and IRAS2A remains cold enough 
for CO to stay frozen, as evidenced by the green region in the 
bottom left panel of Fig. |9] As a result, horizontal cuts from the 
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Fig. 6. Gas temperature (black) and dust temperature (red) along 
the cavity wall of the three sources due to UV heating. Shocks 
are not included. The coordinate on the horizontal axis is the dis- 
tance to the protostar (Eq. (|2]l). The blue curves show the temper- 
ature profiles due to passive heating only (gas and dust coupled). 



in Sect. I4.2l bv fitting a grid of luminosities to the CO observa- 
tions. The resulting temperature profiles along the cavity walls 
are shown in Fig. |6] (black and red curves). For comparison, the 
plot also shows the temperature profiles from DUSTY (blue), 
which is what the gas and dust would reach in the absence of 
UV heating. The effect of UV heating is strongest for DK Cha 
because of the lower envelope densities. The gas temperature 
peaks at 3800 K at 150 AU from the star IRAS2A reaches its 
maximum of 850 K right at the inner edge, followed by a small 
secondary maximum of 330 K at 320 AU. For HH 46, the density 
in the inner envelope is high enough that the gas initially remains 
coupled to the dust and has the same temperature. At larger dis- 
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Fig. 8. Top: radial abundance profile of CO (red) and H2O (blue) 
in the spherical envelope models with passive heating only. 
Bottom: abundance profiles along the cavity wall (Eq. (|2])) in the 
2D models, with all effects of the UV field included: photodis- 
sociation, photodesorption, and altered chemistry due to heating 
of the gas and dust. Shocks are not included in either panel. 



cavity wall into the envelope essentially show drop-abundance 
profiles (Figs. [T2]and[T3]in the online appendix). CO is not ex- 
pected to freeze out in the less massive and warmer envelope 
of DK Cha, so the horizontal cuts merely show the effect of the 
photodissociation rate decreasing with distance to the cavity wall 
(Figs.[T0]and[T4li. 

The protostellar UV field also dissociates H2O, making for 
a nearly uniform abundance of lO '" along the wall of each of 
the three sources (Fig. [8]). Moving into the envelope, the increas- 
ing amount of extinction slows down the photodissociation pro- 
cess and the H2O abundance increases to a few 10 in the two 
regions where it does not freeze out: the warm inner envelope 
and the low-density outer envelope (Fig.|9l bottom right). In be- 
tween, freeze-out keeps the abundance limited to at most a few 
10"^ As in the case of CO, the original ID abundance profiles 
largely survive at the midplane (Figs. fTTHl4l in the online ap- 
pendix). 

With the conditions in the passively heated gas and the UV- 
heated gas established, the final question is what happens in the 
small-scale shocks along the cavity walls. As shown in Fig. [3] the 
peak temperature in a 20 km s"' shock ranges from 1300 to 3800 
K for pre-shock densities from lO'* to lO*"^ cm"^. In addition, the 
peak te mperature depends roughly linearly on the shock veloc- 
ity (e.g. jFlower & Pineau Des Foretsll20"T0l) . The higher temper- 
atures help overcome the barriers on the O -h H2 ^ OH + H 
and OH -H H2 ^ H2O + H reactions (Bergi n et al]ll998h . so 
the H2O abundance in the shocked gas is expected to be orders 
of magnitude higher than the value of 10 '° calculated for the 
UV-heated gas (Fig.[8l). 



02468 10 02468 10 
R (1000 AU) 

Fig. 9. Abundance of CO (left) and H2O (right) in two versions 
of the HH 46 model: the spherical envelope with passive heating 
only (top) and the envelope with UV-heated outflow cavity walls 
(bottom). 

The shock models of i Flower & Pineau des ForetsI (l2003h 
predict a peak H2O abundance of ~10~'* (Fig. |3]l, but that is in 
the absence of a UV field (Sect. l23T l. The photodissociation and 
reformation timescales can be estimated to get an idea of how 
much the H2O abundance would change in an irradiated shock. 
At 1000 AU from the protostar, the UV flux is equivalent to 
Go ~ 100, which gives a photodissociation timescale of about 1 
yr (Ivan Dishoeck et al ][2006). A reformation timescale of about 
1 yr o ccurs for r = 1000 K and n(H2) = 10^ cm"^^ (Bergin et all 
Il998h . The density at 1000 AU is closer to 10^ cm'^ (Fig.©, 
so the reformation timescale is a factor of a few shorter than 
the photodissociation timescale. However, the two numbers are 
similar enough that in reality either process could dominate over 
the other, in particular when considering different points along 
the cavity wall. A new generation of shock models, incorporat- 
ing UV radiation, is needed to solve this issue. Until then, the 
peak H2O abundance of ~10 in Fig.|3]should be considered an 
upper limit only. 

4.2. UV luminosities and sliock velocities 

The model has two free parameters - the UV luminosity and 
the shock velocity - for which the best-fit values are obtained by 
comparing the computed CO line fluxes to the available Herschel 
observations in a procedure. The lower-7 CO lines observed 
from the ground are not included in the fit. They contain contri- 
butions from the parent molecular cloud and from colder swept- 
up gas, neither of which is present in our model. The H2O obser- 
vations are also not used to constrain the model. The grid of UV 
luminosities runs from 0.01 to 1 L^. The grid of shock velocities 
runs from 10 to 35 km s"', approximately the critical velocity of 
C-type shocks for the relevant density regime (Sect. 123] ). 
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Fig. 15. Top: CO fluxes computed for the passive and UV com- 
ponents in HH 46 together, for UV luminosities ranging from 
0.06 to 0. 15 Lq- Bottom: CO fluxes computed for the shock com- 
ponent in HH 46, for shock velocities ranging from 10 to 35 km 



As an example of the dynamic range of the free parameters. 
Fig. [15] shows the predicted line fluxes for HH 46 for a range 
of values. The top panel shows how the combined flux from 
the passively heated envelope and the UV-heated outflow cavity 
walls changes as function of UV luminosity. The bottom panel 
does so for the flux from the small-scale shocks as function of 
shock velocity. These figures are qualitatively the same for the 
other two sources. 

The best-fit UV luminosities and shock velocities are tabu- 
lated in Table |2] and shown in Fig. [16] Overplotted in the figure 
are the 1, 2 and 3cr confidence intervals. Both parameters are 
well constrained in HH 46. In IRAS2A, the UV luminosity is 
only constrained from the top. In DK Cha, the shock velocity is 
unconstrained at the 3<t level and only constrained from below at 
the 2cr level. Extending the grid of luminosities to lower values 
does not resolve the issue for IRAS2A, because the gas cannot 
cool below the dust temperature. The analysis suggests that 
the CO fluxes observed with Herschel are powered primarily by 
shocks in IRAS2A, by UV-heating in DK Cha, and by a combi- 
nation of both in HH 46. This result is discussed in more detail 
in the next section. 

4.3. Line fluxes 

4.3.1. Carbon monoxide 

With the physical and chemical structure from Sect. 14.11 and 
the best-fit UV luminosities and shock velocities from Sect. 
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Fig. 16. Reduced contours for the UV luminosity and the 
shock velocity in each source. Red crosses mark the minima. 
Contours are at the 1,2 and 3cr levels. 



132] (taking Luv = 0.01 Lq for IRAS2A and u, = 35 km s"' 
for DK Cha), the model produces the integrated CO line fluxes 
shown in Fig.[T7] Also plotted are the integrated line fluxes ob- 
served from the ground and with Herschel. The only systematic 
discrepancy between model and observations is an underproduc- 
tion of the 2-1, 3-2 and 4 -3 lines. The same discrepancy was 
seen bv lvan Kempen et"al] ((2009) in their purely spherical model 
of HH 46. The excess flux probably originates in cold gas in the 
parent molecular cloud or in cold swept-up outflow gas, neither 
of which is included in our model. The low-7 lines of the less 
abundant C"^0 isotopologue are not expe cted to have any signif- 
icant cloud or out flow contribution (e.g., |j0rgensen et al.ll2b02l; 
lYildiz et al.ll2010l) . and indeed the model reproduces the avail- 
able C'^O data to within 50%. The mismatches with some of the 
individual Herschel data points in Fig.[T7]do not appear to fol- 
low any systematic trend. They are likely due to uncertainties in 
the model fluxes and minor factors such as pointing off'sets and 
calibration uncertainties that may be larger than the estimated 
10-30%. 

The coloured curves in Fig.[T7]show the predicted contribu- 
tion from the passively heated envelope (green), the UV-heated 
cavity walls (red) and the small-scale shocks (blue) to the to- 
tal CO emission (black). Each source requires a combination of 
UV-heated gas and small-scale shocks to reproduce the Herschel 
data, although there are clear qualitative differences from source 
to source. IRAS2A is the youngest of the three sources (deeply 
embedded Stage 0) and its CO ladder appears to be dominated by 
shock-powered emission. HH 46 is more evolved (Stage I) and 
shows a 50/50 split between UV-powered and shock-powered 
emission. DK Cha is still more evolved (in transition from Stage 
I to II) and its small-scale shocks are responsible for only a small 
fraction of the overall CO emission. The decreasing contribu- 
tion of the shock component is consistent with the general trend 
of the outflow force b eing weaker in more evolved protostars 
(iBontemps et al.|[T996l) . The increasing contribution of the UV 
component is powered by the decreasing densities in the enve- 
lope, allowing the gas to become hotter (Figs.|5]and|6j see also 
Sect. 15. lb . Although the sample size is small, the observed data 
are entirely in keeping with these physically motivated trends. 
Application of the model to a larger number of sources, such as 
the full WISH sample, will help to obtain a firmer conclusion. 
Furthermore, a full parameter study is planned to assess system- 
atically how the CO emission depends on properties such as the 
envelope mass and the UV luminosity. 

The contribution in each curve can be quantified by sum- 
ming the line fluxes from 7u = 2 to 40 and correcting them for 
distance to give the overall CO luminosities for the adopted in- 
clinations (Table O. The CO line fluxes can also be plotted in a 
rotational diagram to investigate excitation conditions. Table [5] 
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Fig. 17. CO line fluxes observed with the JCMT, APEX and Herschel. Overplotted are the fluxes from the three model components 
(green, red and blue) and their sums (black). Error bars (Icr) reflect the 10-30% calibration uncertainty of Herschel and the 20-30% 
calibration uncertainty of the JCMT and APEX. Upper limits are at the 3cr level. 



Table 4. Absolute and relative luminosities of CO and H2O for 
the three physical components in each source. 



Component 


1RAS2A 


HH46 


DK Cha 


for CO 


10-3 ^ 


% 


10-3 Lo 


% 


10-3 Lo 


% 


passive 


0.3 


7 


0.5 


7 


0.1 


2 


UV 


0.8 


20 


3.1 


45 


4.0 


78 


C-shocks 


3.0 


73 


3.3 


48 


1.0 


20 


total 


4.1 


100 


6.9 


100 


5.1 


100 




Component 


1RAS2A 


HH46 


DK Cha 


for H2O 


10-3 


% 


10-3 Lo 


% 


10-3 


% 


passive 


0.02 


0.02 


0.53 


0.38 


0.05 


0.17 


UV 


0.001 


0.001 


0.04 


0.03 


3(-5) 


l(-4) 


C-shocks 


100 


>99.9 


140 


99.6 


30 


99.8 


total 


100 


100 


140 


100 


30 


100 



Notes. The tabulated values are for the adopted source inclinations of 
90, 53 and 0° (Table |2ll. They do not represent the total line cooling 
budget, which would require integrating over all inclinations from edge- 
on to pole-on. 



Table 5. Rotational temperatures for CO and H2O in the three 
physical components in each source. 



Component 


1RAS2A 


HH 


46 


DK 


Cha 


for CO 




(K) 




(K) 




(K) 


passive 


31 


±5 


33 


±9 


54 


±5 


UV 


240 


±30 


160 


±20 


410 H 


= 100 


C-shocks 


380 


±20 


400 


±20 


1600 


±50 




Component 


1RAS2A 


HH 


46 


DK 


Cha 


for H2O 




(K) 




(K) 




(K) 


passive 


72 


± 1 


130 


±5 


175 


± 10 


UV 


30 


± 1 


69 


± 1 


114 


± 1 


C-shocks 


410 


±80 


540 H 


t 170 


750 H 


= 400 



lists the rotational temperatures derived by fitting straight lines 
to the passive curve from - 10 to 100 K (up to 7 = 5-4), to 
the UV curve from 100 to 1000 K (up to 7 = 18-17), and to the 
shock curve from 1000 to 4000 K. The error margins reflect the 
stochastic Icr uncertainty on each fit, but do not account for any 
uncertainties in the model results. 



The CO ladder for DK Cha in Fig.[niimphes that the six CO 
lines seen with ISO, from J - 14-13 to 19-18 (iGiannini et al.1 
[T99l . originate in the UV-heated gas in the outflow cavity 
walls. Based on a pass ively heated spherical envelope model, 
Ivan Kempen et al.l (f2006) concluded that the ISO lines trace ma- 
terial on larger spatial scales than the material responsible for 
the 7-6 and lower lines observed with APEX. With the more 
complex source geometry in our model, such a requirement no 
longer exists. In fact, the temperature gradient along the cavity 
wall (Fig.|6]l, together with the pole-on inclination, suggests the 
high-7 lines seen with ISO originate on smaller spatial scales 
than the lower-/ APEX lines. A more detailed analysis of the 
spatial extent of the PACS data, using all 5 x 5 spaxels, can help 
to solve the puzzle, but this is beyond the scope of the current 
work. 

4.3.2. Water 

Figure [18] shows the observations and model predictions for 
H2O. The complicated rotational structure of H2O does not al- 
low for plots in the format of the CO ladders in Fig. [17] so the 
results are shown as rotational diagrams instead. The individual 
line fluxes can be summed again to obtain the total line luminosi- 
ties. However, the models only calculate fluxes for a subset of all 
possible rotational lines. Fluxes for the missing lines in the HIFI 
and PACS wavelength ranges (up io E^lk - 2000 K) are cal- 
culated by fitting a single rotational temperature to the available 
lines in each model component (Table O. The level populations 
of H2O are not in LTE, so this method may produce large errors 
in the flux of a given individual line. However, because those 
errors tend to cancel when summing over all possible lines, as- 
suming a single temperature i s appropriate for calculatin g the 
total line luminosities (see also [Goldsmith & Langeill999l) . The 
resulting values are tabulated in Table [4] 

Unlike CO, the H2O emission observed with Herschel seems 
to be dominated (>99%) by shock-heating in all three sources. 
This is consistent with the broad H2O line p rofiles seen with 
HIFI in IRAS2A and other embedded YSOs dKristensen et al.l 
12010 '. subm.). Our shock models actually overproduce many of 
the detections and upper limits by up to a factor of 10, in partic- 
ular in HH 46 (Fig. [TST l. The likely cause is that the H2O abun- 
dances in our shock models were calculated without taking the 
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Fig. 18. Rotational diagrams for water Black: Herschel detections (circles) and 3cr upper limits (triangles). Green, red and blue: 
model predictions for the passive, UV and shock components. 



protostellar UV field into account. Based on the photodissocia- 
tion and reformation timescales estimated in Sect. 14.11 the H2O 
abundance can easily be lowered by a factor of 2 to 10 from the 
peak value of lO"'* shown in Fig.|3] At the same time, the abun- 
dance of OH would go up, so our shock models are currently ex- 
pected to underproduce the OH line intensities. Ho wever, as ar - 
gued by .van Kempen et al.. (2010b ) and W ampfler et al.l (1201 Oh . 
the OH lines observed with Herschel Ukely originate in a differ- 
ent physical component than the H2O lines (a dissociative J-type 
shock instead of a non-dissociative C-type shock). It is therefore 
not possible to test whether the OH line intensities predicted for 
the small-scale shocks are indeed lower than the actual OH emis- 
sion from that component. A potential problem with lowering 
the H2O abundances through photodissociation is that some CO 
in the shock may be dissociated as well. If that happens, a higher 
shock velocity is required to fit the CO observations, which in 
turn would raise the model H2O fluxes back to their original val- 
ues of up to an order of magnitude too high. An independent con- 
straint on the CO and H2O abundances in UV-irradiated shocks 
is needed to fully resolve this issue. 

Another difference between CO and H2O is seen in the rela- 
tive contributions of the passive and UV components. For H2O, 
both are weak compared to the shock-powered emission, but the 
passive component is always the brighter of the two (Table |4|. 
The weak emission from the UV-heated gas in the outflow cav- 
ity walls is due to the very low H2O abundance (Figs.[8]and|9]), 
which in turn is due to the wide range of photon energies over 
which H2O can be dissociated. CO can only be dissociated by 
photons of at least 1 1 eV, so its abundance in the cavity walls re- 
mains higher and the UV-heated gas contributes a larger fraction 
of the overall CO emission. 



5. Discussion 

5.1. Gas temperature 

The largest unknown in our model is the gas temperature due 
to UV heating of the outflow cavity walls. Calculating the gas 
and dust temperatures for a given density and UV flux is a stan- 
dard problem in models of photon-dominated regions (PDRs; 
e.g., Tielens & Hollenbach 1985; Hollenbach & fielens 1997). 
However, the PDR community has not yet converged on a 
unique solution: in a benchmark test of ten different PDR codes. 



the temperature solutions varied by up to an order of magni- 
tude (Rolligetal. 2007) . As shown in Sect. 14.31 the tempera- 
ture grid from Kaufma n et alJ (1 19991 hereafter K99), along with 
the exp(-0.6Av) depth dependence derived from the benchmark 
test, works well to reproduce the CO observations. However, 
given the complicated nature of our model, this does not neces- 
sarily mean that the,K99 grid gives the correct temperatures. For 
example, if the grid systematically overpredicts the temperature, 
it may still give a good match to the data by decreasing the UV 
luminosity. Other temperature grids may also be able to match 
the data. The goal of this section is two-fold: (1) to show that, 
at least on a qualitative level, the .K99 grid follows the density 
and UV flux dependence expected from basic heating and cool- 
ing physics for the relevant range of conditions; and (2) to see 
how robust our results are to changes in the adopted temperature 
grid. 

The EH grid is reproduced in Fig. \T9\ In its original form, 
it only covers H2 densities up to 10^ cm"^. The densities in 
our model go up to 10"^ cm"-', so the grid is extrapolated as 
Tgas l/n(H2). The n"' dependence comes from the approxi- 
mate functional forms of the heating and cooling mechanisms 
expected to dominate at high densities: photoelectric heating 
(oc n; Bakes & Tielens 1994) and gas-grain collisional cooling 
(oc rp-; IHollenbach & M cKee 1989). The density dependence of 
the temperature in the original .K99. grid is almost exactly n ' 
between 10^ and 10^ cm"-'. 

Overplotted as black lines in Fig.[T9]are the gas density and 
incident UV flux along the cavity walls of the three sources. The 
density and UV flux both fall off approximately as r"^, so they 
scale linearly with each other At densities above ~10^ cm"^, the 
gas temperature is expected to scale with the i nverse of the den- 
sity ( see above) and linearly with the UV flux (iBakes & TielensI 
Il994i) . and should therefore remain constant along the inner part 
of the wall. This is indeed seen in Fig. |6] Tgas does not change 
much over the range of radii corresponding to densities of more 
than 10^ cm"^. Photoelectric heating likely remains dominant at 
lower densities, but atomic fine-structure lines take over from 
gas-grain collisions as the dominant cooling mechanism. Their 
cooling rate depends linearly on the density as long as the den- 
sity exceeds the critical density of 3 x lO-' cm""* for the [Cii] 
158 fim line, so the te mperature is expected to re main roughly 
constant with density (IK99I ; iMeijerink etani2007h . Because the 
Go dependence remains, the temperature along the cavity wall 
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Fig. 19. PDR surface temperatu re grid (colou r scale and solid 
white contours) from Kaufman et alJ (Il999h . extrapolated as 
Tgas 0^ n ' for n(H2) > 10^ cm"^. The white contours are drawn 
at 10, 30, 100 and 1000 K. Overplotted as dashed white contours 
(same values) is the approximate temperature grid based on Eqs. 
^ and (|7]). The black lines trace the conditions along the cavity 
walls of the three sources. 
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Fig. 20. CO ladder for HH 46. The data points are as in Fig. [H] 
The dashed black curve is the model flux for the U V compo- 
nent for temperatures from the 'Kaufman e t alj (1 19991) grid. The 
dashed coloured curves are the same, but for the temperature 
approximation from Eqs. (|6]l and (|7ll, using different heating ef- 
ficiencies. The solid red curve combines the dashed red curve 
with the passive and shock curves from Fig. [17] 



should decrease as r ' for n(H2) < 10^ cm"^. Again, this is the 
approximate behaviour seen in Fig.|6l The dependencies change 
again for densities below 3 x lO"* cm"-', but such low densities 
do not appear in our model. 

The location of the black lines in Fig.[T9lshows why the UV- 
powered CO emission looks qualitatively different in DK Cha 
than in IRAS2A and HH 46. The lower envelope densities in the 
former turn the conditions along the cavity wall into a classic 
PDR, with densities of 10''-10^ cm^^ and UV fluxes of 10^- 
10"* times the mean interstellar flux. The gas therefore gets much 
hotter and CO is excited to higher rotational levels. The steep 
temperature gradient in this part of Fig. [19] also means that the 
CO emission is rather sensitive to details in the SED fitting with 
DUSTY. A factor of two difference in density leads to a simi- 
lar change in the gas temperature, which can affect the highest- 
] CO lines by a few orders of magnitude. Figure [T7] tentatively 
shows an evolutionary trend of the CO emission in more evolved 
sources being more strongly UV-powered. While the application 
of a physical model is necessary to quantify exactly how much 
of the total CO ladder in a given source is due to UV heating, 
the observations are consistent with the main point that with 
decreasing density (assuming sufficient UV), sources becomes 
more FDR-like. This is what one would naturally expect in the 
transition from Stage to Stages I and II objects as presented in 
this paper 

Of equal importance to the CO emission is the choice of the 
PDR temperature grid. An easy way to get a rough idea of how 
changes in the adopted grid would affect the CO line fluxes is 
to parameterise the gas temperature based on the dominant heat- 
ing and coohng mechanisms. Two equations suffice to do so. 
The first one approximates th e gas-grain collisional cooling rate 
(iHollenbach & McKed[T989l) : 



A,_, = (lxlO-ergcm-s-)4y^^- 
75 K\^ 



,-3 „-K„2 



T 



100 A 



1 - 0.8 exp 



T 



(7"ga.s ~ 7"(]ust) . 



(6) 



where nn is the total hydrogen nucleus density (equal to 2n(H2) 
in our model for all practical purposes) and amin is the minimum 
grain size (set to 10 A). The second equation is for the photo- 
electric heating rate dBakes & Tielensll 19941) : 



- pe 



(1 X 10"^'*ergcm"^ s"')eGonH ■ 



(7) 



The photoelectric heating efficiency, e, is treated as a free pa- 
rameter for the purpose of this simple approximation. As such, 
it hides any uncertainties in the grain charge an d in the nu- 
merical prefactors in the two equations (see, e.g., lYoung et al.l 
2004). Given a dust temperature from RA DMC (or, alternat ively, 
from the analytical expression of Hollen bach et aljri99lh . the 
gas temperature at niW-i) > 10^ cm"-' is obtained by solving 
for Ag_g = Fpe. At lower densities, the dependence of Tgas on 
n(H2) disappears to zeroth order, and we set rgas[Go, n(H2) < 
10*^ cm^] = Tg^,[GQ,lQ^ cm^]. 

Figure [19] shows the 10, 30, 100 and 1000 K contours for 
e = 0.2 a s dashed white lines. The approximate grid reproduces 
the 'K9S' grid to within a factor of 2 for the full range of densi- 
ties and UV fluxes encountered along the three cavity walls. The 
approximation breaks down for Go/n(H2) < lO"'* (high density, 
low UV) and Go/n(H2) > 10"^ (low density, high UV), and can- 
not be used in either o f thos e regimes. The value of 0.2 required 
for a good match with K99 exceeds the maximum value of 0.05 
allowed by Bakes & Tielens ( 1994). However, the approximate 
temperature grid is purely meant as an easy tool to study the ef- 
fect of different temperature grids, so this factor of 4 shoul d not 
be considered physically important. The key point is that the lK99l 
grid and the approximate grid show the same qualitat ive tr ends 
for the relevant range of conditions, which means the lK99l grid 
follows the behaviour expected from the dominant heating and 
cooling mechanisms. 

With £ - 0.2, the approximate grid does a good job of re- 
producing the CO line emission for the UV component. Shown 
in Fig. [20] are the integrated line fluxes from the UV-heated gas 
in HH 46, using either the K99 grid (dashed black curve, identi- 
cal to the red curve from Fig. [TTl i or the approximate grid with 
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three different values of e (d ashed coloured curves). The e = 0.2 
curve lies within 30% of the curve up to 7u = 20, at which 
point shocks take over as the dominant producer of CO emis- 
sion. When adding the passive and shock components to the 
e = 0.2 curve (producing the solid red curve), the overall agree- 
ment with the data is as good as in Fig. [T7] Taking a lower or 
higher photoelectric heating efficiency would under- or overpro- 
duce the J = 14-13, 16-15 and 18-17 lines. 

The photoelectric heating rate scales linearly with e and Go 
(Eq. (|7|i), so the three values of e in Fig. |20]can be brought to 
give the same CO fluxes by changing the UV luminosity. Herein 
lies the degeneracy alluded to at the beginning of this section 
and in Sect. 14.21 as long as two temperature grids show the same 
qualitative trends, both can reproduce the CO observations, but 
they have different best-fit UV luminosities. This degeneracy 
can only be lifted by calibrating the temperature grids against 
a source whose incident UV flux can be constrained indepen- 
dently in the same part of n(H2)-Go parameter space. 



5.2. Caveats in the shock models 

The shock modelling relies on three assumptions: (1) the shocks 
move along the cavity walls, (2) the effective beam-filling fac- 
tor can be obtained from accurately estimating the shock width, 
and (3) the protostellar UV field does not affect the H2O abun- 
dance. In order to test where the emission is coming from in this 
scenario - extended lower-density gas with a large beam-filling 
factor or higher-density gas with a small beam-filling factor - 
fluxes from the standard shock model are compared to fluxes 
from models without any contribution from gas with densities of 
10"* cm"^ and lower, 10^ cm"-' and lower, or lO*" cm"-' and lower. 
The shock velocity is kept at 20 km s"' for this experiment. For 
the case of HH 46, gas at densities below 10"* cm"-' is found to 
contribute less than 10""*%, and gas at densities between 10"* and 
10^ cm"^ only 0.3%. About 75% of the overall flux comes from 
the segment with a pre-shock density of 10^^ cm"^. The other 
two sources show almost the same numbers, indicating that the 
shocked gas observed by PACS is dominated by emission from 
the densest parts of the YSO. 

If the density is assumed to be fixed at more than 10^ cm"^, 
the shock velocity alone determines the absolute line fluxes. This 
was already seen in the bottom panel of Fig. [15] where the basic 
shape of the CO ladder is nearly independent of velocity. Hence, 
the velocity determines the intensity of the CO emission, and the 
pre-shock density determines the shape of the ladder. Comparing 
the fluxes amongst the three sources, the differences are domi- 
nated by their respective distances. This confirms that the actual 
geometry only plays a small role with respect to the predicted 
CO emission, and that the bulk of emission is created in very 
dense gas contained within the central PACS spaxel. 

If the effective beam-filling factor is treated as a free param- 
eter, it becomes degenerate with the shock velocity. A smaller 
beam-filling factor would require a higher shock velocity to pro- 
duce the same line fluxes, and vice versa. However, for a given 
source geometry, such as the ellipsoid outflow cavity in our 
model, the beam-filling factor follows from the cooling lengths 
and is not truly a free parameter. The only way to constrain ei- 
ther the geometry or the beam-filling factor would be to have 
an independent handle on the pre-shock density and the shock 
velocity. Our method therefore gives a unique solution for the 
given shock geometry, but other geometries cannot be excluded. 



6. Conclusions 

The Herschel Space Observatory offers a new probe into the hot 
gas present in embedded low-mass young stellar objects (YSOs) 
through emission lines of highly rotationally excited CO (up to 
4000 K above the ground state) and H2O (up to 600 K). This 
paper presents the first model to reproduce these observations 
quantitatively. The model starts with a passively heated spheri- 
cal envelope, which fits the lowest-7 CO isotopologue lines ob- 
served from the ground, but falls short of the ^^CO Herschel 
data by several orders of magnitude. Two additional model com- 
ponents are therefore invoked: (1) gas in the walls of a bipo- 
lar outflow cavity, exposed to the protostellar UV field, and (2) 
small-scale C-type shocks along the cavity walls. 

The main conclusion is that, within the context of the model, 
the far-infrared CO rotational line emission in the low-mass 
YSOs NGC1333 IRAS2A, HH 46 and DK Cha can only be 
reproduced with a combination of UV-heated gas and C-type 
shocks. The three sources show a tentative evolutionary trend: 
the CO emission appears to be dominated by shock-heating in 
the youngest source (IRAS2A), by UV-heating in the oldest 
source (DK Cha), and by a 50/50 combination in the source of 
intermediate age (HH 46). This trend is powered by the envelope 
density decreasing as a protostar evolves from Stage to Stages 
I and II. 

The H2O lines are dominated by shocks in all three sources, 
with less than 1 % of the total H2O line luminosity originating in 
UV- and passively heated gas. This is a direct result of the dif- 
ferent abundances in the three components. Freeze-out keeps the 
H2O abundance below 10"** in the cold outer envelope. The dust 
grains in the UV-heated layer are warm enough for H2O to evap- 
orate, but it is dissociated by the UV field once it does. Efficient 
reformation at high temperatures boosts the abundance to 10"^- 
lO"'* in the C-type shocks, so they are responsible for most of 
the H2O emission. CO has a lower evaporation temperature (20 
versus 100 K), so its abundance is orders of magnitude higher 
than that of H2O in the passively heated and the UV-heated gas, 
allowing these two components to contribute more strongly to 
the overall CO emission. 

The main uncertainty in the model lies in the calculation of 
the gas temperature as function of the incident UV flux. This is a 
well-known problem in photon-dominated regions (PDRs), and 
many solutions are available in the literature. However, they dis- 
agree by up to an order of magnitude and some do not even agree 
qualitatively on how the gas temperature depends on the density 
or the UV flux. The sensitivity of our model results to the treat- 
ment of the gas temperature is tested with an approximate ana- 
lytical description of the expected dominant heating and cooling 
mechanisms. This exercise shows that the quantitative contribu- 
tions of the UV-heated gas and the C-type shocks to the hot gas 
emission observed by Herschel are uncertain by at least a fac- 
tor of 2, but the qualitative conclusion that both components are 
needed to explain the observations is robust. 
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NGC1333 IRAS2A HH 46 DK Cha 




R (100 AU) 

Fig. 10. Abundance profiles for CO for the model with outflow cavities and UV heating. Shocks are not included. The bottom row 
zooms in on the regions indicated in the top row. 



NGC1333 IRAS2A HH 46 DK Cha 




R (100 AU) 



Fig. 11. As Fig. [To] but for H2O. 
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Fig. 12. Gas temperature (red) and abundances of CO and H2O (solid and dotted black) in eight horizontal cuts through the envelope 
ofNGC1333 IRAS2A. 
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Fig. 13. Gas temperature (red) and abundances of CO and H2O (solid and dotted black) in eight horizontal cuts through the envelope 
ofHH 46. 
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Fig, 14. Gas temperature (red) and abundances of CO and H2O (soUd and dotted black) in eight horizontal cuts through the envelope 
ofDKCha. 



